1 Transects

Plant Flowers Date lon lat ele Year Month
Glossoloma oblongicalyx 4 2015-10-19 -78.59093 0.130838 2270 2015 10
Gasteranthus quitensis 2 2016-10-17 -78.59770 0.120070 1940 2016 10
Kohleria affinis 1 2016-12-13 -78.59534 0.126746 2110 2016 12
Columnea ciliata 3 2014-02-27 -78.59934 0.116682 1960 2014 2
Columnea medicinalis 1 2014-04-23 -78.59372 0.128700 2130 2014 4
Drymonia teuscheri 3 2016-07-28 -78.59245 0.129393 2200 2016 7

2 Interactions

2.1 Species elevation ranges

2.2 Flowering Data Matrix

3 Are flowering patterns non-random with respect hummingbird overlap and elevation?

Let’s start with a traditional randomization test to get everyone comfortable with the essential result. For each random draw, calculate the mean niche overlap in hummingbird usage. Species can only bloom at site which they occur. Mantains plant abundance.

3.1 Prepare data for JAGS

4 Occurrence

Species elevation models without conditional probability of flowering

## sink("model/occurrence.jags")
## cat("
##     model {
## 
##     for (x in 1:Nobs){
## 
##     #Observation of a flowering plant
##     Y[x] ~ dbern(p[x])
##     logit(p[x]) <-  alpha[Plant[x]] + beta[Plant[x]] * elevation[x]
## 
##     #Residuals
##     discrepancy[x] <- abs(Y[x] - p[x])
## 
##     #Assess Model Fit
##     Ynew[x] ~ dbern(p[x])
##     discrepancy.new[x]<-abs(Ynew[x] - p[x])
##     }
## 
##     #Sum discrepancy
##     fit <- sum(discrepancy)/Nobs
##     fitnew <- sum(discrepancy.new)/Nobs
## 
##     #Prediction
## 
##     for(x in 1:Npreds){
##     #predict value
## 
##     #Observation - probability of flowering
##     prediction[x] ~ dbern(p_new[x])
##     logit(p_new[x]) <- alpha[NewPlant[x]] + beta[NewPlant[x]] * new_elevation[x]
## 
##     #predictive error
##     pred_error[x] <- abs(Ypred[x] - p_new[x])
##     }
## 
##     #Predictive Error
##     fitpred<-sum(pred_error)/Npreds
## 
##     #Priors
## 
##     #Species level priors
## 
##     for (j in 1:Plants){
##       #Intercept flowering probability
##       alpha[j] ~ dnorm(0,0.386)
##       beta[j] ~ dnorm(0,0.386)
##     }
## 
##     }
##     ",fill=TRUE)
## 
## sink()
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 3600
##    Unobserved stochastic nodes: 4416
##    Total graph size: 46411
## 
## Initializing model

5 Get Chains

5.0.1 Evaluate convergence

5.0.2 Posterior estimates

5.0.3 Predict flowering rates by elevation

5.1 One example data point

Confirm everything matches up

## # A tibble: 1 x 8
##   Plant             Site Month Year      n   ele Index jPlant
##   <fct>            <dbl> <dbl> <chr> <dbl> <dbl> <chr>  <dbl>
## 1 Columnea ciliata     3     3 2018      1  1.78 50         2
## # A tibble: 4 x 7
##   Plant             Site Month Year      n   ele Index
##   <fct>            <dbl> <dbl> <chr> <dbl> <dbl> <int>
## 1 Columnea ciliata     3     3 2014      1  1.77    18
## 2 Columnea ciliata     3     3 2015      1  1.76    66
## 3 Columnea ciliata     3     3 2018      1  1.78   201
## 4 Columnea ciliata     3     3 2019      1  1.79   247
##   Draw chain              par     value   parameter jPlant
## 1    1     1 discrepancy[100] 0.6163511 discrepancy      4
## 2    2     1 discrepancy[100] 0.7356973 discrepancy      4
## 3    3     1 discrepancy[100] 0.6173596 discrepancy      4
## 4    4     1 discrepancy[100] 0.6357632 discrepancy      4
## 5    5     1 discrepancy[100] 0.5818248 discrepancy      4
## 6    6     1 discrepancy[100] 0.6252359 discrepancy      4
##                  Plant Index Site Month Year n   ele
## 1 Columnea mastersonii   100    1    12 2016 1 1.364
## 2 Columnea mastersonii   100    1    12 2016 1 1.364
## 3 Columnea mastersonii   100    1    12 2016 1 1.364
## 4 Columnea mastersonii   100    1    12 2016 1 1.364
## 5 Columnea mastersonii   100    1    12 2016 1 1.364
## 6 Columnea mastersonii   100    1    12 2016 1 1.364

Probability of observing a flowering record

Discrepancy

##        mean     upper     lower
## 1 0.6507852 0.7559506 0.5431877

Fit by presence/abs

6 Occurrence + Flowering

Probability of flowering conditional on elevation occurrence

## sink("model/occurrence_month.jags")
## cat("
##     model {
## 
##     for (x in 1:Nobs){
## 
##     #Observation of a flowering plant
##     Y[x] ~ dbern(p[x])
##     logit(p[x]) <-  alpha[Plant[x],Month[x]] + beta[Plant[x]] * elevation[x]
## 
##     #Residuals
##     discrepancy[x] <- abs(Y[x] - p[x])
## 
##     #Assess Model Fit
##     Ynew[x] ~ dbern(p[x])
##     discrepancy.new[x]<-abs(Ynew[x] - p[x])
##     }
## 
##     #Sum discrepancy
##     fit <- sum(discrepancy)/Nobs
##     fitnew <- sum(discrepancy.new)/Nobs
## 
##     #Prediction
## 
##     for(x in 1:Npreds){
##     #predict value
## 
##     #Observation - probability of flowering
##     prediction[x] ~ dbern(p_new[x])
##     logit(p_new[x]) <- alpha[NewPlant[x],NewMonth[x]] + beta[NewPlant[x]] * new_elevation[x]
## 
##     #predictive error
##     pred_error[x] <- abs(Ypred[x] - p_new[x])
##     }
## 
##     #Predictive Error
##     fitpred<-sum(pred_error)/Npreds
## 
##     #Priors
##     #Species level priors
##     for (j in 1:Plants){
##     #effect of elevation
##       beta[j] ~ dnorm(0,0.386)
##     for(k in 1:Months){
##       #Intercept flowering probability by month
##       alpha[j,k] ~ dnorm(0,0.386)
##     }
## 
##     }
## 
##     }
##     ",fill=TRUE)
## 
## sink()
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 3600
##    Unobserved stochastic nodes: 4592
##    Total graph size: 51196
## 
## Initializing model

7 Get Chains

7.0.1 Evaluate convergence

Elevation response

month intercepts

7.0.2 Posterior estimates

Fit by presence/absence

7.1 One example data point

Confirm everything matches up

## # A tibble: 1 x 8
##   Plant                Site Month Year      n   ele Index jPlant
##   <fct>               <dbl> <dbl> <chr> <dbl> <dbl> <chr>  <dbl>
## 1 Besleria solanoides     5     9 2014      1  2.14 20         1
## # A tibble: 5 x 7
##   Plant                Site Month Year      n   ele Index
##   <fct>               <dbl> <dbl> <chr> <dbl> <dbl> <int>
## 1 Besleria solanoides     5     9 2014      1  2.14    41
## 2 Besleria solanoides     5     9 2015      1  2.14    88
## 3 Besleria solanoides     5     9 2016      1  2.14   126
## 4 Besleria solanoides     5     9 2017      0  2.23   173
## 5 Besleria solanoides     5     9 2018      0  2.23   225

Probability of observing a flowering record

Discrepancy

##       mean     upper     lower
## 1 0.684529 0.8883212 0.4516712

7.1.1 Predict flowering rates

7.1.2 Predict elevation ranges

8 Predictive Models of Co-Flowering

8.1 Visitor Attraction

9 Get Chains

9.1 Attraction

## sink("model/occurrence_attraction.jags")
## cat("
##     model {
##     
##     for (x in 1:Nobs){
##     #Observation of a flowering plant
##     Y[x] ~ dbern(p[x])
##     logit(p[x]) <-   e[Plant[x],Month[x]] + beta[Plant[x]] * elevation[x]
##     
##     #Residuals
##     discrepancy[x] <- abs(Y[x] - p[x])
##     
##     #Assess Model Fit
##     Ynew[x] ~ dbern(p[x])
##     discrepancy.new[x]<-abs(Ynew[x] - p[x])
##     }
##     
##     #Sum discrepancy
##     fit <- sum(discrepancy)/Nobs
##     fitnew <- sum(discrepancy.new)/Nobs
##     
##     #Prediction
##     
##     for(x in 1:Npreds){
##     #predict value
##     
##     #Observation - probability of flowering
##     prediction[x] ~ dbern(p_new[x])
##     logit(p_new[x]) <- e[NewPlant[x],NewMonth[x]] + beta[NewPlant[x]] * new_elevation[x] 
##     
##     #predictive error
##     pred_error[x] <- abs(Ypred[x] - p_new[x])
##     }
##     
##     #Predictive Error
##     fitpred<-sum(pred_error)/Npreds
##     
##      #########################
##     #autocorrelation in error
##     #########################
##     
##     #For each of month, error covariance among flowering plants
##     for(k in 1:Months){
##       e[1:Plants,k] ~ dmnorm(zeros,tauC[,])
##     }
## 
##     ##covariance among similar species
##     for(i in 1:Plants){
##     for(j in 1:Plants){
##     C[i,j] = exp(-lambda_cov* D[i,j])
##     }
##     }
##     
##     ## Covert variance to precision for each parameter, allow omega to shrink to identity matrix
##     vCov = omega*C[,] + (1-omega) * I
##     tauC = inverse(vCov*gamma^2)
##     
##     #Autocorrelation priors
##     gamma  ~ dunif(0,20)
##     
##     #Strength of covariance decay
##     lambda_cov ~ dunif(0.1,2)
##     omega ~ dbeta(1,1)
## 
##     #Priors
##     #Species level priors
##     for (j in 1:Plants){
##       #effect of elevation
##       beta[j] ~ dnorm(0,0.386)
##       alpha[j] ~ dnorm(0,0.386)
##     }
##     
##     }
##     ",fill=TRUE)
## 
## sink()
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 3600
##    Unobserved stochastic nodes: 4431
##    Total graph size: 52012
## 
## Initializing model

9.1.1 Evaluate convergence

Elevation response

month intercepts

9.1.2 Posterior estimates

9.2 Effect of interaction attraction

9.3 Repulsion

## sink("model/occurrence_repulsion.jags")
## cat("
##     model {
##     
##     for (x in 1:Nobs){
##     #Observation of a flowering plant
##     Y[x] ~ dbern(p[x])
##     logit(p[x]) <-  e[Plant[x],Month[x]] + beta[Plant[x]] * elevation[x]
##     
##     #Residuals
##     discrepancy[x] <- abs(Y[x] - p[x])
##     
##     #Assess Model Fit
##     Ynew[x] ~ dbern(p[x])
##     discrepancy.new[x]<-abs(Ynew[x] - p[x])
##     }
##     
##     #Sum discrepancy
##     fit <- sum(discrepancy)/Nobs
##     fitnew <- sum(discrepancy.new)/Nobs
##     
##     #Prediction
##     
##     for(x in 1:Npreds){
##     #predict value
##     
##     #Observation - probability of flowering
##     prediction[x] ~ dbern(p_new[x])
##     logit(p_new[x]) <-  e[NewPlant[x],NewMonth[x]] + beta[NewPlant[x]] * new_elevation[x]
##     
##     #predictive error
##     pred_error[x] <- abs(Ypred[x] - p_new[x])
##     }
##     
##     #Predictive Error
##     fitpred<-sum(pred_error)/Npreds
##     
##     #########################
##     #autocorrelation in error
##     #########################
##     
##     #Error covariance among flowering plants
##     for(k in 1:Months){
##       e[1:Plants,k] ~ dmnorm(zeros,tauC[,])
##     }
##     
##     ##covariance among similar species
##     for(i in 1:Plants){
##     for(j in 1:Plants){
##     C[i,j] = exp(-lambda_cov*D[i,j])
##     }
##     }
##     
##     ## Covert variance to precision for each parameter, allow omega to shrink to identity matrix
##     vCov = omega*C[,] + (1-omega) * I
##     tauC = vCov*gamma^2
##     
##     #Autocorrelation priors
##     gamma  ~ dunif(0,20)
##     
##     #Strength of covariance decay
##     lambda_cov ~ dunif(0.1,2)
##     omega ~ dbeta(1,1)
##     
##     #Priors
##     #Species level priors
##     for (j in 1:Plants){
##       #effect of elevation
##       beta[j] ~ dnorm(0,0.386)
##       alpha[j] ~ dnorm(0,0.386)
##     }
##     
##     }
##     ",fill=TRUE)
## 
## sink()
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 3600
##    Unobserved stochastic nodes: 4431
##    Total graph size: 52011
## 
## Initializing model

10 Get Chains

10.0.1 Evaluate convergence

10.0.2 Posterior estimates

10.1 Effect of interaction repulsion

11 Model Comparison

11.0.1 e: The effect of autocorrelation

Glossoloma purpureum

11.1 Alpha: Species specific flowering rate

11.2 Omega: The magnitude of the effect of autocorrelation on mean flowering occurrence

11.3 Gamma: The variance of the effect of autocorrelation on mean flowering occurrence

11.4 Lambda: The decay in autocorrelation effect

11.5 Decay in autocorrelation effect

12 Model Fit

12.1 Bayesian pvalue

## # A tibble: 4 x 2
##   Model                      p
##   <chr>                  <dbl>
## 1 baseline               0.854
## 2 interaction_attraction 0.599
## 3 interaction_repulsion  0.485
## 4 occurrence_month       0.42

Without baseline 

## # A tibble: 3 x 2
##   Model                      p
##   <chr>                  <dbl>
## 1 interaction_attraction 0.599
## 2 interaction_repulsion  0.485
## 3 occurrence_month       0.42

12.2 Model Fit

Model mean lower upper
interaction_repulsion 0.1894753 0.1839611 0.1954302
occurrence_month 0.1846309 0.1789345 0.1904198
interaction_attraction 0.1685083 0.1628142 0.1742645

12.2.1 Without baseline

12.3 By Species

12.3.1 Without baseline

12.3.2 Zoom in

13 Prediction

13.1 Example from Glossoloma oblongicalyx

13.1.1 Tables

Model mean lower upper
baseline 0.1847199 0.1783802 0.1911261
interaction_repulsion 0.1769840 0.1705929 0.1838919
occurrence_month 0.1734432 0.1665216 0.1802577
interaction_attraction 0.1660734 0.1592176 0.1728697

14 Flowering correlations in predictions

15 Example

Columnea medicinialis v Columnea strigosa which have strong overlap in visitors

Dint["Columnea medicinalis","Columnea strigosa"]
## [1] 0.3083944

Logit e: Autocorrelation among species

InvLogit e: Autocorrelation among species

## # A tibble: 2 x 4
## # Groups:   Model [2]
##   Model                  Var1                Var2             Correlation_E
##   <chr>                  <fct>               <fct>                    <dbl>
## 1 interaction_attraction Columnea medicinal… Columnea strigo…            NA
## 2 interaction_repulsion  Columnea medicinal… Columnea strigo…            NA
## # A tibble: 4 x 4
## # Groups:   Model [4]
##   Model                  Var1                Var2             Correlation_P
##   <chr>                  <fct>               <fct>                    <dbl>
## 1 baseline               Columnea medicinal… Columnea strigo…         0.991
## 2 interaction_attraction Columnea medicinal… Columnea strigo…         0.963
## 3 interaction_repulsion  Columnea medicinal… Columnea strigo…         0.963
## 4 occurrence_month       Columnea medicinal… Columnea strigo…         0.774
## # A tibble: 4 x 4
## # Groups:   Model [4]
##   Model                 Var1               Var2            Correlation_Ynew
##   <chr>                 <fct>              <fct>                      <dbl>
## 1 baseline              Columnea medicina… Columnea strig…           NA    
## 2 interaction_attracti… Columnea medicina… Columnea strig…            0.768
## 3 interaction_repulsion Columnea medicina… Columnea strig…           -0.136
## 4 occurrence_month      Columnea medicina… Columnea strig…           -0.630

16 Figures for the paper

17 Figure 1

Flowering patterns

18 Example Flowering patterns

19 Figure 4

Model Fit and Prediction Fit

20 Predictions from the baseline models

Comparison of the attraction and repulsion

20.1 For a given plant

20.1.1 Observed

Overlay observed

Prediction

Predictions from the best model

Supplamental figures